#####################################
#Stewart & Zhukov (09)
#Civil Military Article Figures and
#Descriptive Statistics on Multiply Imputed DV
#THIS CODE IS LARGELY UNCOMMENTED BUT SHOULD RUN CLEANLY GIVEN FILE SET
#####################################
rm(list=ls(all=TRUE))
library(car)
library(MASS)
set.seed(12345)
##Data Merging and TimeSeries Operations
probmatrix <- read.delim('ProbabilityMatrix.csv', sep=',', header=T)
topics <- read.delim('TopicModels.csv', sep=',', header=T)
data <- read.delim('TFIDFLSA.csv', sep=',')
probmatrix <- cbind(data$labels, probmatrix)
rm(data)
m1 <- merge(topics, probmatrix, by.x = "filename", by.y = "data$labels")
OtherData <- read.delim('090206_cm.csv', sep=',', header=T) 
m2 <- merge(m1, OtherData, by.x = "DOCUMENT_ID", by.y = "DOCUMENT_ID", all.x=FALSE, all.y=FALSE)  
rm(OtherData)
rm(topics)

head(probmatrix)
problist <- cbind(probmatrix$wp1, probmatrix$wp2, probmatrix$wp3)
draws <- function(probs) {
  sample(x=c(1,2,3), size= 10000, prob=probs, replace=TRUE)	
}
multiclass <- apply(problist, 1, draws)  
multiclass <- t(multiclass)
rm(draws)
rm(probmatrix)

########
#Figure 1
########
  means <- function(column) sum(column==1)
  samplingdist <- apply(multiclass, 2, means)
  quantile(samplingdist, probs=c(.05, .95)); mean(samplingdist)
#4166, 4208.753, 4251

  means <- function(column) sum(column==2)
  samplingdist <- apply(multiclass, 2, means)
  quantile(samplingdist, probs=c(.05, .95)); mean(samplingdist)
#886, 915.7715, 945 

  means <- function(column) sum(column==3)
  samplingdist <- apply(multiclass, 2, means)
  quantile(samplingdist, probs=c(.05, .95)); mean(samplingdist)
#2754, 2791.476, 2830 




m1$topics <- recode(m1$topics, "0 = 'NA'; 1='NA'; 2=1; 3='NA'; 4=0; 5=0; 6=0")
m2$topics <- recode(m2$topics, "0 = 'NA'; 1='NA'; 2=1; 3='NA'; 4=0; 5=0; 6=0")
m2$topics <- as.numeric(m2$topics)
sum(na.omit(m1$topics==1))
sum(na.omit(m1$topics==0))
4291

classbardata <- c(4208.753, 915.7715, 2791.476)
topicbardata <- c(1923, 1702, 4291)

par(mfrow=c(1,2))

barplot2(topicbardata, names.arg=c("Interventionist","Realpolitik","None of the Above"),width=c(1),space=0,ylim=c(0,5000),
col="lightgrey",ylab="# of Documents",main="Issue Salience by Document", density=c(10))

barplot2(classbardata, names.arg=c("Conservative","Activist","None of the Above"),width=c(1),space=0,ylim=c(0,5000),
col="lightgrey",ylab="# of Documents",main="Use of Force by Document", density=c(10))
lower <- c(4166, 886, 2754)
upper <- c(4251,945,2830)
segments(x0=c(1:3)-.5, y0=lower, x1=c(1:3)-.5, y1=upper,lwd=10, col="black", lend="butt")
legend(x="topright",lty=c("solid"),
legend=c("95% Confidence Intervals"),lwd=c(3),col=c("black"))

##########
#Figure 2
##########
par(mfrow=c(1,1))
##Upper Left
attach(m2)
cornerUL1 <- sum(na.omit(topics[MILT==1 & INFORMAL==0]))
cornerUL1[2] <- sum(na.omit(topics[MILT==0 & INFORMAL==0]))
cornerUL1[3] <- sum(na.omit(topics[MILT==1 & INFORMAL==1]))
cornerUL1[4] <- sum(na.omit(topics[MILT==0 & INFORMAL==1]))
cornerUL1[5] <- sum(na.omit(topics[PMNT==1]))
cornerUL1[6] <- sum(na.omit(topics[CABT==1]))
cornerUL2 <- sum(na.omit(topics[MILT==1 & INFORMAL==0]==0))
cornerUL2[2] <- sum(na.omit(topics[MILT==0 & INFORMAL==0]==0))
cornerUL2[3] <- sum(na.omit(topics[MILT==1 & INFORMAL==1]==0))
cornerUL2[4] <- sum(na.omit(topics[MILT==0 & INFORMAL==1]==0))
cornerUL2[5] <- sum(na.omit(topics[PMNT==1]==0))
cornerUL2[6] <- sum(na.omit(topics[CABT==1]==0))
cornerUL <- cbind(cornerUL1, cornerUL2)

mybarcol <- "gray20"
mp <- barplot2(t(cornerUL), beside = TRUE,
        col = c("darkgrey", "lightgrey"),
        legend = c("Interventionist Statements", "Realpolitik Statements"), ylim = c(0, 1500),
        main = "Issue Salience by Group Statements (Total Number of Documents)", font.main = 40, names.arg=c("Military (formal)","Political (formal)","Military (informal)",
"Political (informal)","Parliament","Cabinet"), ylab="Number of Statements",
        plot.grid = TRUE)
mtext(side = 1, at = colMeans(mp), line = 2,
      text = paste("Mean", formatC(rowMeans(cornerUL))), col = "black")
box()


##Upper Right
cornerUR1 <- sum(na.omit(topics[MILT==1 & INFORMAL==0]))/sum(MILT==1&INFORMAL==0)
cornerUR1[2] <- sum(na.omit(topics[MILT==0 & INFORMAL==0]))/sum(MILT==0 & INFORMAL==0)
cornerUR1[3] <- sum(na.omit(topics[MILT==1 & INFORMAL==1]))/sum(MILT==1 & INFORMAL==1)
cornerUR1[4] <- sum(na.omit(topics[MILT==0 & INFORMAL==1]))/sum(MILT==0 & INFORMAL==1)
cornerUR1[5] <- sum(na.omit(topics[PMNT==1]))/sum(PMNT==1)
cornerUR1[6] <- sum(na.omit(topics[CABT==1]))/sum(CABT==1)
cornerUR2 <- sum(na.omit(topics[MILT==1 & INFORMAL==0]==0))/sum(MILT==1&INFORMAL==0)
cornerUR2[2] <- sum(na.omit(topics[MILT==0 & INFORMAL==0]==0))/sum(MILT==0 & INFORMAL==0)
cornerUR2[3] <- sum(na.omit(topics[MILT==1 & INFORMAL==1]==0))/sum(MILT==1 & INFORMAL==1)
cornerUR2[4] <- sum(na.omit(topics[MILT==0 & INFORMAL==1]==0))/sum(MILT==0 & INFORMAL==1)
cornerUR2[5] <- sum(na.omit(topics[PMNT==1]==0))/sum(PMNT==1)
cornerUR2[6] <- sum(na.omit(topics[CABT==1]==0))/sum(CABT==1)
cornerUR <- cbind(cornerUR1, cornerUR2)
detach(m2)
mybarcol <- "gray20"
mp <- barplot2(t(cornerUR), beside = TRUE,
        col = c("darkgrey", "lightgrey"),
        legend = c("Interventionist Statements", "Realpolitik Statements"), ylim = c(0, 1), ylab="Proportion of Group Statements",
        main = "Issue Salience (Proportion by Group)", font.main = 40, names.arg=c("Military (formal)","Political (formal)","Military (informal)",
"Political (informal)","Parliament","Cabinet"), sub="Number of Documents (Excluding NAs)",
        plot.grid = TRUE)
mtext(side = 1, at = colMeans(mp), line = 2,
      text = paste("N=", formatC(rowSums(cornerUL))), col = "black")
box()

cornerUR

#############
##Bottom Left
#############
attach(m2)
  means <- function(column) sum(column==1 & test==1)

###BL1
test <- c(MILT==1 & INFORMAL==0)
  samplingdist <- apply(multiclass, 2, means)
ci1.1 <- quantile(samplingdist, probs=c(.05))
ci1.2 <- quantile(samplingdist, probs=c(.95))
cornerBL1 <- mean(samplingdist)

test <- c(MILT==0 & INFORMAL==0)
  samplingdist <- apply(multiclass, 2, means)
ci1.1[2] <- quantile(samplingdist, probs=c(.05))
ci1.2[2] <- quantile(samplingdist, probs=c(.95))
cornerBL1[2] <- mean(samplingdist)

test <- c(MILT==1 & INFORMAL==1)
  samplingdist <- apply(multiclass, 2, means)
ci1.1[3] <- quantile(samplingdist, probs=c(.05))
ci1.2[3] <- quantile(samplingdist, probs=c(.95))
cornerBL1[3] <- mean(samplingdist)

test <- c(MILT==0 & INFORMAL==1)
  samplingdist <- apply(multiclass, 2, means)
ci1.1[4] <- quantile(samplingdist, probs=c(.05))
ci1.2[4] <- quantile(samplingdist, probs=c(.95))
cornerBL1[4] <- mean(samplingdist)

test <- c(PMNT==1)
  samplingdist <- apply(multiclass, 2, means)
ci1.1[5] <- quantile(samplingdist, probs=c(.05))
ci1.2[5] <- quantile(samplingdist, probs=c(.95))
cornerBL1[5] <- mean(samplingdist)

test <- c(CABT==1)
  samplingdist <- apply(multiclass, 2, means)
ci1.1[6] <- quantile(samplingdist, probs=c(.05))
ci1.2[6] <- quantile(samplingdist, probs=c(.95))
cornerBL1[6] <- mean(samplingdist)

###BL2
  means <- function(column) sum(column==2 & test==1)

test <- c(MILT==1 & INFORMAL==0)
  samplingdist <- apply(multiclass, 2, means)
ci2.1 <- quantile(samplingdist, probs=c(.05))
ci2.2 <- quantile(samplingdist, probs=c(.95))
cornerBL2 <- mean(samplingdist)

test <- c(MILT==0 & INFORMAL==0)
  samplingdist <- apply(multiclass, 2, means)
ci2.1[2] <- quantile(samplingdist, probs=c(.05))
ci2.2[2] <- quantile(samplingdist, probs=c(.95))
cornerBL2[2] <- mean(samplingdist)

test <- c(MILT==1 & INFORMAL==1)
  samplingdist <- apply(multiclass, 2, means)
ci2.1[3] <- quantile(samplingdist, probs=c(.05))
ci2.2[3] <- quantile(samplingdist, probs=c(.95))
cornerBL2[3] <- mean(samplingdist)

test <- c(MILT==0 & INFORMAL==1)
  samplingdist <- apply(multiclass, 2, means)
ci2.1[4] <- quantile(samplingdist, probs=c(.05))
ci2.2[4] <- quantile(samplingdist, probs=c(.95))
cornerBL2[4] <- mean(samplingdist)

test <- c(PMNT==1)
  samplingdist <- apply(multiclass, 2, means)
ci2.1[5] <- quantile(samplingdist, probs=c(.05))
ci2.2[5] <- quantile(samplingdist, probs=c(.95))
cornerBL2[5] <- mean(samplingdist)

test <- c(CABT==1)
  samplingdist <- apply(multiclass, 2, means)
ci2.1[6] <- quantile(samplingdist, probs=c(.05))
ci2.2[6] <- quantile(samplingdist, probs=c(.95))
cornerBL2[6] <- mean(samplingdist)

cornerBL <- cbind(cornerBL1, cornerBL2)
ci.1 <- cbind(ci1.1, ci2.1)
ci.2 <- cbind(ci1.2, ci2.2)

normweights <- matrix(c(sum(MILT==1&INFORMAL==0),sum(MILT==0 & INFORMAL==0),sum(MILT==1 & INFORMAL==1),sum(MILT==0 & INFORMAL==1)
  ,sum(PMNT==1),sum(CABT==1)), nrow=6, ncol=2)
detach(m2)
cornerBR <- cornerBL/normweights
ciw.1 <- ci.1/normweights
ciw.2 <- ci.2/normweights

par(mfrow=c(1,1))

mybarcol <- "gray20"
mp <- barplot2(t(cornerBL), beside = TRUE,
        col = c("darkgrey", "lightgrey"),
        ylim = c(0, 3000),
        main = "Attitudes Towards Use of Force by Group (Total Number of Documents)", font.main = 40, sub = "95 percent error bars", col.sub = mybarcol,
        plot.ci = TRUE, ci.l = t(ci.1), ci.u = t(ci.2), names.arg=c("Military (formal)","Political (formal)","Military (informal)",
"Political (informal)","Parliament","Cabinet"), ylab="Number of Statements",
        plot.grid = TRUE)
mtext(side = 1, at = colMeans(mp), line = 2,
      text = paste("Mean", formatC(rowMeans(cornerBL))), col = "black")
box()

mybarcol <- "gray20"
mp <- barplot2(t(cornerBR), beside = TRUE,
        col = c("darkgrey", "lightgrey"),
        legend = c("Conservative Statements", "Activist Statements"), ylim = c(0, 1),
        main = "Attitudes Towards Use of Force (Proportion by Group)", font.main = 40, sub = "95 percent error bars", col.sub = mybarcol,
        plot.ci = TRUE, ci.l = t(ciw.1), ci.u = t(ciw.2), names.arg=c("Military (formal)","Political (formal)","Military (informal)",
        "Political (informal)","Parliament","Cabinet"), ylab="Proportion of Group Statements",
        plot.grid = TRUE)
mtext(side = 1, at = colMeans(mp), line = 2,
      text = paste("Avg. N", formatC(rowSums(cornerBL))), col = "black")
box()







##############################
#Figure 3
#SPATIAL PLOT
##############################

dates <- unique(m2$YRMO[m2$topics==1|m2$topics==0])
dates <- dates[-16]




naomitmean <- function(item) mean(na.omit(item))
yprime <- aggregate(m2$topics, list(Date = m2$YRMO, Military = m2$MILT), FUN=naomitmean) 
yprime <- yprime[c(-5,-9, -34, -128),]

xP <-c(0)
yP <-c(0)
xM <-c(0)
yM <-c(0)

for (i in 1:3) {
  z <-  recode(multiclass[,i], "1=0; 2=1; 3=NA")
  z <- aggregate(z, list(Date = m2$YRMO, Military = m2$MILT), FUN=naomitmean)
  z <- merge(z, yprime, by.x = c("Date", "Military"), by.y = c("Date", "Military"))
  xP <- append(xP, z[z$Military==0,3])
  yP <- append(yP, z[z$Military==0,4])
  xM <- append(xM, z[z$Military==1,3])
  yM <- append(yM, z[z$Military==1,4])
}

plot(xP, yP, pch=22, ylab="0=Realpolitik, 1=Interventionist", xlab="0=Conservative, 1=Activist", 
  main="Group Disagreement: 3 Simulated Time-Series Runs", cex.lab=1.3, cex.main=1.3)
points(xM, yM, pch=19)
legend(x="top",pch=c(22,19),
legend=c("Political", "Military"))


#############################
##Figure 4
##Disagreement
#############################

xprime <- aggregate(m2$SC.MILP, list(Date = m2$YRMO), FUN=mean) 
pprime <- aggregate(m2$PFSCR, list(Date = m2$YRMO), FUN=mean)
ix <- c()
dy <- c()

sims <- 20
coefmat <- matrix(NA, nrow=sims, ncol=3)
for (i in 1:sims) {
  z <-  recode(multiclass[,i], "1=0; 2=1; 3=NA")
  z <- aggregate(z, list(Date = m2$YRMO, Military = m2$MILT), FUN=naomitmean)
  z <- merge(z, yprime, by.x = c("Date", "Military"), by.y = c("Date", "Military"))
  z <- merge(z, xprime, by.x=c("Date"), by.y=c("Date"))
  z <- merge(z, pprime, by.x=c("Date"), by.y=c("Date"))
  z <- na.omit(z)
  mil <- z[z$Military==1,]
  pol <- z[z$Military==0,]
  z <- merge(mil, pol, by.x=c("Date", "x.x.1","x.y.1"), by.y=c("Date", "x.x.1", "x.y.1"))
  z <- na.omit(z)
  disagreement <- sqrt((z$x.x.x - z$x.x.y)^2 + (z$x.y.x - z$x.y.y)^2)
  coefmat[i,] <- coef(lm(disagreement ~ z$x.x.1 + z$x.y.1))
  ix <- append(ix, z$x.x.1)
  dy <- append(dy, disagreement)
}

summary(lm(disagreement ~ z$x.x.1 + z$x.y.1))

plot(jitter(ix, factor=2), dy, main="Elite Integration over 20 Simulations", xlab="Proportion of Military Elites on Security Council", 
  ylab="Disagreement (Spatial Metric)", 
  xlim=c(0,max(ixnew)), sub="Controlling for Press Freedom Score")
for (i in 1:sims) {
  abline(a=(coefmat[i,1] + coefmat[i,3]*mean(PFSCR)), b= coefmat[i,2], col="grey")
}













































